Variable-range projection model for turbulence-driven collisions 
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We discuss the probability distribution of relative speed AV of inertial particles suspended in a 
highly turbulent gas when the Stokes numbers, a dimensionless measure of their inertia, is large. We 
identify a mechanism giving rise to the distribution P(AV) ~ exp(— C| AV| 4//3 ) (for some constant 
C). Our conclusions are supported by numerical simulations and the analytical solution of a model 
equation of motion. The results determine the rate of collisions between suspended particles. They 
are relevant to the hypothesised mechanism for formation of planets by aggregation of dust particles 
in circumstellar nebula. 

PACS numbers: 05.20.Dd, 45.50.Tn, 47. 27.-i, 47.57.E- 
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1. Introduction. It is widely believed that the first 
stage of the formation of planets involves the aggregation 
of microscopic dust grains in the gaseous nebula around 
young stars This process must occur in a turbulent 
environment, because the transport of angular momen- 
tum by diffusion would be too slow to account for the 
lifetimes of these nebula. Also, the aggregation process 
occurs in gas with a very low density, so that the motion 
of the dust grains is very lightly damped. It is necessary 
to achieve a good understanding of the relative veloc- 
ity of collisions of the dust grains to determine whether 
and how planet formation could result from the aggre- 
gation of microscopic dust grains. The relative velocity 
is required to determine the rate of collision of the dust 
grains. Also, if the relative velocity is sufficiently high, 
clusters may fragment upon collision. These issues con- 
cerning planet formation are discussed in 0, Q • 

Earlier discussions of the relative velocity of suspended 
particles @, [E @| have estimated the order of magnitude 
of the relative velocity, but a satisfactory theory for its 
distribution has been lacking. In the context of planet 
formation, the case of lightly damped particles is most 
important. If the microscopic correlation time of the flow 
is r and the damping rate (defined by (|2|) below) is 7, we 
define the Stokes number as St = I/7T. A theoretical ap- 
proach is required, because simulations are impracticable 
for the lightly damped case where St ;§> 1. 

In this letter we show that the probability distribu- 
tion function for the relative velocities AV of colliding 
particles is well approximated by 



P(AV) = Aexp f-C|AF| 4 / 3 7 2/3 /£ 2/3 



(1) 



where Z is the turbulence intensity (the rate of dissipation 
per unit mass) and C is a universal dimensionless constant 
(with A determined by normalising the distribution). We 
argue that this is a precise asymptote for the distribution 
for large AT^|. 




FIG. 1: Variable-range projection model. We show curves of 
constant probability p(Ax,Av) for Ax S> 1 (black), At; oc 
Ax 1 / 3 . Also shown is a realisation of a trajectory of equation 
© projected from large separations to Aa; = (red), com- 
pared with optimal trajectories (blue dashed). Here Ax and 
Av are dimensionless variables in equation (|SJ. 



We remark that there are connections with the distri- 
bution of accelerations in turbulent flows. The accelera- 
tion of a suspended particle is proportional to its velocity 
relative to the fluid. Because the relative velocity of two 
particles with St 3> 1 is the sum of their (statistically in- 
dependent) velocities relative to the fluid, the tail of the 
distribution of accelerations a of suspended particles is 
of the form P(a) ~ exp[— const|a| 4 / 3 ], analogous to |T]). 
For suspended particles with St -c 1, the acceleration is 
the same as Langrangian fluid acceleration, which also 
has a distribution of the same form as Jl}, with 4/3 re- 
placed by « 2/5 0- The distribution of accelerations for 
suspended particles in a turbulent flow was studied nu- 
merically for a range of values of St by Bee et al Q . The 
results (figure 2b of their paper) are compatible with the 
limiting cases discussed above. 

Our explanation of the mechanism underlying equation 
(HJ) proceeds as follows. The colliding particles acquire 
a relative velocity when they are accelerated by differ- 
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FIG. 2: Probability density p(Ax, Av) for a simulation of 
equations ©, (o), compared with theory (solid line): a 
Ax = 0, compared with fTJ. b Ax = 2.5, compared with 
distribution obtained from (|9I10I11[) . In both cases a = 2/3 
and e — 1. The empirical distributions are normalised and 
the prefactor of the theoretical curves chosen to give the best 
fit. Panels c and d are the same as for a and b respectively, 
except a = 4/3. 



ent regions of the fluid. They are then 'projected' (i.e. 
thrown) a certain distance away from the fluid element 
which accelerated them. Since relative particle veloci- 
ties imparted by the fluid flow increase with separation, 
particles which collide with a high relative velocity ac- 
quired their relative motion when their separation was 
large. Our estimate of the probability distribution func- 
tion P(AV) involves a maximisation of the probability 
of reaching zero separation with respect to variation of 
the distance over which the particles are projected by the 
flow. We term this model the 'variable range projection' 
model. It has much in common with the 'variable-range 
hopping' model for electrical conduction in semiconduc- 
tors at low temperatures Q, which also arises from an 
optimisation of the hopping length and leads to an ex- 
pression for the conductance of the form (fTJ) , with tem- 
perature playing the role of the relative velocity. 

Our heuristic description is supported by precise 
asymptotic analysis of a one-dimensional model, equa- 
tion ^ below. Figure [5h shows a comparison with sim- 
ulation. For non-zero separation, the relative velocity 
distribution has a more complex asymmetric form, Fig. 
[2]b. We also confirm a surmise about the variance of the 
relative velocity @ . 

2. Equations of motion. The equations of motion for 
the position x and velocity v of a suspended particle are 
r = v and v = j[u(r,t) — v] , where u(r,t) is the fluid 
velocity. This equation is applicable even when the gas 
mean free path is large compared to the size of the par- 



ticles [10(. The corresponding equation for the relative 
displacement AX and relative velocity AV of two par- 
ticles is 

AX = AV , AV = 7 [Am(AX, t) - AV] (2) 

and where Am = u(AX,t) — u(0,t). According to the 
Kolmogorov theory of turbulence, there is a range of 
lengthscales £ for which a component Au of the relative 
velocity of fluid elements with separation £ is determined 
only by the turbulence intensity. Dimensional arguments 
[111 ] then imply 

(Au(£, t)Au{£, 0)} = («) 2 / 3 /(t£ 1/3 /^ 2/3 ) (3) 

for some function / (angular brackets are used to denote 
averages throughout this paper). 

3. Variable-range projection model. Consider the rel- 
ative displacement AX and speed AV of two particles. 
When AX is small, the driving effect of the fluid velocity 
Au is negligible, and the damping term is most signifi- 
cant. At greater distances, the relative velocity of the 
background fluid drives the relative motion of the parti- 
cles. First let us consider the relative motion in greater 
detail at small separations, such that we can neglect the 
effect of the driving term Am. In this case AV decays ex- 
ponentially in time, so that if two particles collide with 
relative velocity AV at time t, their relative velocity at 
an earlier time to was AVo(to) = AV exp[y(t — to))]. In- 
tegrating this expression, we find that the relative sepa- 
ration at time to was 



AAo(io) 



dt' AVe^'^ — d 

1 



1 -7(*-*o) 



(4) 

so that AVo — AV — 7AA0, where AAo was the ini- 
tial separation. Continuing to neglect the effects of the 
fluid velocity, we see that in order for particles to collide 
with relative velocity AV, they must have had a larger 
velocity difference AVo a t a larger, and unknown, separa- 
tion AX . For large AX , equations resemble those of 
an Ornstein-Uhlenbeck process [jjj , where the velocity is 
Gaussian distributed. We therefore expect that for suffi- 
ciently large AXq, the relative velocity is approximately 
Gaussian distributed: 



p(AV ,AX ) 



1 



exp 



2<Al/ 2 ) 



(5) 



Here we use the expectation that for large separations, 
the relative velocity is well approximated by the relative 
velocity of the fluid elements, so that equation ([3]) im- 
plies that (AV^ 2 ) ~ (£ AX ) 2/3 . To determine where the 
inbound particle colliding with relative velocity AV orig- 
inated, we therefore find the value of the separation AAo 
which maximises the probability of colliding with relative 
velocity AV, that is we maximise p(AVo, AXq), where 
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AV = AV-jAX , with respect to AX . Figure [T] illus- 
trates the trajectories. Let the value for which the max- 
imum obtains be AXq. Neglecting the pre-exponential 
factor of ©, we find AA * = -AV/2j. The distribu- 
tion of velocities for colliding particles is predicted to be 
P(AV) = p(AV - j AXq, AX*). Neglecting the pre- 
exponential factor, we obtain equation ([I]). For the vari- 
ance of the relative velocity, it follows that (Ay 2 ) oc £/7- 
This provides a justification for a result which was previ- 
ously inferred from the Kolmogorov theory of turbulence 
by a dimensional argument [6(. 

4. Microscopic model. The motion of the smallest ed- 
dies in a fully-developed turbulent flow are characterised 
by the Kolmogorov length 77, Kolmogorov time r and Kol- 
mogorov velocity i*k- The flow is characterised by two 
dimensionless variables, the Stokes number, St = I/7T, 
and the Kubo number, Ku = ukt/tj. In turbulent veloc- 
ity fields 77, t and uk are functions of the dissipation rate 
£ and the kinematic viscosity v. Dimensional consider- 
ations then imply that Ku = 0(1). However in the fol- 
lowing we consider a model for a turbulent flow in which 
Ku <C 1 , corresponding to a very rapidly fluctuating flow 
field, which can be modelled by a Langevin equation. 

Consider the equations of motion (|2|) in one spatial di- 
mension. We convert to dimensionless variables, writing 
t' = jt, Ax = AX/rj, Av — AV/rj"/. When the velocity 
field Au is very rapidly fluctuating, we can approximate 
the equation of motion in scaled variables by the follow- 
ing Langevin equation 

dAx = Avdt', dAv = -Avdt' + Sw (6) 

where the random increment Sw satisfies 

(8w) =0, (Sw 2 ) = 2V(Ax)dt' , V{Ax) = e|Ax| Q . (7) 

Here we have introduced a parameter e ~ Ku 2 . Hav- 
ing approximated equations © by equations ^ , we find 
that solutions of ((6]) for different values of e can be ob- 
tained from the solution with e = 1 by a scaling transfor- 
mation. Although it suffices to consider the case where 
e = 1, we retain e in subsequent expressions because it 
will be used as a small parameter of a WKB expansion. 
This formal procedure allows us to study the tails of the 
joint probability distribution of Ax and Av in a con- 
trolled manner. In (|7|) we also allow for an arbitrary 
exponent < a < 2. The value of a is determined by re- 
quiring that the variance of the relative velocity has the 
correct behaviour as Ax — > 00: the solution of Q pre- 
sented below indicates that (Av 2 ) ~ |Ax| Q for Ax — > 00, 
so comparison with (J3j) indicates that a = 2/3 is the 
correct choice. 

5. Distribution of collision velocities. The distribu- 
tion |1| of collision velocities is determined by the joint 
distribution p(Ax,Av) evaluated at Ax = 0. To deter- 
mine p(Ax, Av) we solve the steady-state Fokker-Planck 
equation corresponding to equations ©, ((7|): 



We note that at large values of Ax (Ax 3> Au), the 
distribution p(Ax,Av) is Gaussian in Av [see equation 
JE}]. In order to solve ^Sj) we make a WKB ansatz [3] 

p(Ax, Av) = K(Ax, Av) exp[-S , (Ax, Av)/e] . (9) 



We write 



S(Ax,Av) = \Ax\ 2 - a g (z, Ax) 
K(Ax,Av) = exp[-gi(z,Ax)], 



(10) 
(11) 



where z — s\Av/Ax (s± = ±1 is chosen so that z > 
0). Assuming that g~o(Ax,z) does not depend on Ax, 
substituting (fT0)l . (fTTj) into |[8]), and collecting terms in 
e" 1 , we obtain 



z(si + z) + s 2 Jz 2 (z + si) 2 - 4g (z)zsi(2 - a) 

ffo(z) = 2^ 

(12) 

where S2 = ±1 labels which branch of the square root is 
to be chosen. In the following we label the solutions of 
(p~2|) by g^ 1 (z). Which of the solutions must be picked 
is determined by the boundary conditions. 

Let us first consider an initial condition (Aar, Av) with 
a positive and large value of Ax. Since z > by def- 
inition, si determines the sign of Av. At large values 
of Ax we know that the distribution of Av is Gaussian 
[eq. (J5J)]. This determines the small- z asymptote of go : 
S = Atj 2 /(2|Ax| q ) = |Ax| 2 ~ Q z 2 /2. Thus we must re- 
quire go ~ 2 2 /2 as z — > 0. We find that only the 
solutions #q ' ' and g^ + match this boundary con- 
dition. In order to reach Aa; = from Aa; > the 
initial relative velocity must be negative. For Ax > 
we are thus forced to choose si = —1, that is to con- 
sider the branch '~ . Consider the case depicted 
in Fig. Q] of a particle projected to Ax = deter- 
mining the distribution of collision velocities. The ac- 
tion is determined by the large- z behaviour of g$ ' \ 
that is S = limAx^o \Ax\ 2 - a g Q ^^ ) (-Av/Ax). We 
find gfi~'~\z) ~ ao( a ) z 2 ~ a for large z. The prefactor 
ao(a) is determined by numerical integration. We find 
do(2/3) w 0.870. The resulting action at Ax = is 



S(Ax = 0, Aw) = a (a)Av 2 



(13) 



To determine the prefactor consider terms of order e° 
arising from substituting (fT0|) . (fTTj) into p]): 

= g' ' — 1 - s\xzd x g\ + (z + s\z 2 - 2g' ) d z g\ (14) 

We make the following separation ansatz g\(x, z) — 
A log Ax + gi(z). It is motivated by the fact that it al- 
lows us to match p(Ax, Aw) to the known behaviour ([5]) 
at large separations. Inserting this ansatz into (fT4]) we 
obtain (neglecting a normalisation constant) 



= - Av Oax P + 9av (Av p) 



e\Ax\ a d 2 Av p 



(8) 



gi = A log Ax 



dz ,l- 9 ^) + Sl Xz' 



z' + siz' 
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Consider now the limiting form of the prefactor K for 
large and for small separations Ax. First, the limit of 
large Ax corresponds to the limit z — > 0. In this limit 
gi is constant and to match the prefactor to the known 
behaviour ([5]) we must set A = 3a/2. Second, the limit 
of Ax — > corresponds to the limit of z — » oo. In this 
limit the integrand in (|15p behaves as ~ X/z' = 3a/ (2z r ). 
Integrating over z we find that e - ^ 1 = |Av|~ 3q / 2 . The 
final result (neglecting a normalisation factor) is thus 

p(0, Av) = |Aw|- 3Q / 2 exp [ - e^aoia^Av] 2 - ] . (16) 

This result, for a — 2/3, corresponds to the distribu- 
tion ([I]) predicted by the variable-range projection model. 
But here it has been derived, including the algebraic pref- 
actor, from a microscopic model. Fig. [2] a, c compares 
of (fl6|) with simulations of the Langevin equation ([6]). 

6. Relative velocities at larger separations. For non- 
zero separations, our WKB approximation is complicated 
by the fact that different branches, corresponding to dif- 
ferent choices signs s±, S2 in (|12p , must be combined. 
For each branch, at finite values of Ax, the contribu- 
tion to p(Ax, Av) is of the form ([5]), with the action 
given by (fit))) and l[T2"|) , with the prefactor given by eqs. 
(fTTjl and (|15p . Which branches must be chosen depends 
upon the signs of Ax and Av. If two branches con- 
tribute for given values of Ax and Aw, the branch with 
the smallest action dominates. The branches which are 
available correspond to four different choices of signs in 
the construction of solutions of (fT2| . namely gQ Sl ' S2 \z). 
We already noted that only the solutions g$ ' \z) and 

.9o + \ z ) can match the correct asymptotic behaviour at 
small z, namely go ~ z 2 /2. 

Let us consider the case where Ax > 0. When Av < 
(that is, when s\ = —1), we find that only the branch 
with action determined by the function ' (z) con- 
tributes, with corresponding action 

S(Ax, Av) = lAx^g^^i-Av/Ax) . (17) 

This expression tends to (TIB")) as Av — > — oo, and to the 
Gaussian form S(Ax, Av) ~ e _1 Av 2 /\ Ax\ a for small val- 
ues of Av. 

For Av > however, the WKB solution is more com- 
plicated. For small z, and for sufficiently small Av 
the solution is given by the branch g { +:+ \z). This 
solution increases very rapidly as z increases; we find 
9o ' ( z ) ~ (1 + a)z 3 /9 as z — * oo, so this branch of 
the WKB solution becomes very small for large Av. By 
adapting the argument in section 3 above, however, we 
can argue that the tails of the probability density for the 
velocity should in fact be given by a branch where the 
action is S ~ ao(a)Av 2 ~ a for Av — > oo, where the pref- 
actor do (a) is the same as for the Av < branch. It is 
possible to find a solution for the branch <7q + ' \z) with 
the correct behaviour, namely <7q + ' \z) ~ ao(a) z 2 ~ a as 



z —> oo. This condition also ensures that the tails of 
p(Ax, Av) are consistent with (fl3|) in the limit Ax — > 0. 

For Av > 0, we therefore construct the solution using 
two branches. For < z < z* the solution constructed 
from go' + \z), satisfying the (z) ~ z 2 /2 for z — > 0, 

is dominant. For z > z*, the solution constructed from 
gQ + '~\z), satisfying g^ + '~\z) ~ a (a)z 2_Q dominates. 
The point z* is determined by the condition that the 
action of the two solutions is equal, that is g^ + ' + \z*) — 
gg + ' \z*). We remark that the solution '~\z) only 
exists for z > z c , where z c is the critical point at which 
the discriminant in (|12[) vanishes. Fortunately, we find 
z* > z c (for a — 2/3 we find z c pa 0.14). The prefactor is 
given by eqs. (jTTJ) and (TT5]). Figure d compares our 
distribution with simulations for Ax ^ 0. 

7. Conclusions. In this letter we have shown how the 
distribution of relative velocities of particles suspended 
in highly turbulent flow at large St = 1/(jt) may be 
surmised from an optimisation argument which we term 
'variable range projection', leading to equation ([1]). We 
validated this simple and general heuristic argument by 
a WKB analysis of a one-dimensional Langevin equation 
model, which produces an identical relative velocity dis- 
tribution at zero separation. 
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